Elastic behavior in Contact Dynamics of rigid particles 
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The systematic errors due to the practical implementation of the Contact Dynamics method for 
simulation of dense granular media are examined. It is shown that, using the usual iterative solver 
to simulate a chain of rigid particles, effective elasticity and sound propagation with a finite velocity 
occur. The characteristics of these phenomena are investigated analytically and numerically in order 
to assess the limits of applicability of this simulation method and to compare it with soft particle 
molecular dynamics. 

PACS numbers: 02.70.Ns,45.70.Cc 



I. INTRODUCTION 

In computational physics one can distinguish two dif- 
ferent validation tasks, which have to be solved in order 
to make simulations a useful research tool. First one must 
prove the validity of a simulation model by comparing its 
results to laboratory experiments, and second, equally 
important, one must assess the systematic errors due to 
the practical implementation in order to tell, how pre- 
cisely the simulation results reflect the theoretical prop- 
erties of the simulation model. In this paper we address 
the second type of validation problem for the simula- 
tion technique of contact dynamics, which was developed 
about 10 years ago 0, ||, || with the aim to investigate 
granular media M in the limit of high rigidity of the par- 
ticles at a high packing density. This method has been 
successfully applied and reproduces experiments (see e.g. 
[^). However, its systematic errors and the computa- 
tional effort to keep them tolerably small have not been 
investigated in detail before. 

This is in marked contrast to other discrete elements 
methods for granular media (cf. e.g. |^), in particular 
the soft particle molecular dynamics simulation model, 
which has been widely used since more than 20 years . 
During this time possible pitfalls such as the detachment 
effect^ and the brake failure effect^ could be discov- 
ered, analyzed and hence avoided. 

Contact dynamics simulations have been applied to 
study a large variety of questions in dense granular sys- 
tems, where excluded volume interactions and static fric- 
tion, so called unilateral constraints, are believed to be 
essential @ ||, |l|| , but it should be mentioned 
that such constraints arise also in other areas like vir- 
tual reality, engineering, especially in robotics, and oper- 
ations research, where the numerical treatments are sim- 
ilar n 0. 

In a system of perfectly rigid particles the sound veloc- 
ity would be infinite. This is in principle borne out by the 
contact dynamics simulation model. However, its practi- 
cal implementation will normally give rise to sound-waves 
in the granular material, as we are going to show in the 



following, even if each single collision is modeled as being 
perfectly inelastic. Our aim is to elucidate this artifact 
after introducing briefly the principles of this method. 



II. THE CONTACT DYNAMICS METHOD 

First, let us point out the basic difference between 
molecular dynamics (MD) on one side and contact dy- 
namics (CD) on the other. Both have in common the 
integration of Newton's equation of motion where the oc- 
curring forces are due to external fields (gravity) or are - 
more important - contact forces, i.e. caused by contacts 
between particles or their contacts with confining walls. 

The spirit of MD is to calculate the contact forces 
according to their cause, i.e. the (usually microscopic) 
deformation of the contact region and the involved ve- 
locities. Since the full treatment of every particle as a 
deformable body would render the simulation of a large 
number of particles exceedingly time consuming, a lot of 
models exist, how to replace the deformation by the local 
overlap [T^ , the latter being a virtual quantity obtained 
from the undeformed shapes. 

The principles of CD are different: Here the contact 
forces are calculated by virtue of their effect, which is 
to fulfill certain constraints. Typically, such a constraint 
is the volume exclusion of the particles or the absence of 
sliding due to static friction. As can be seen immediately, 
this problem cannot be solved locally: In a cluster of par- 
ticles where many contacts are simultaneously present, 
the force at one contact depends on adjacent contacts 
and so on. In that case the aim is to find a global force- 
network, which is consistent with the constraints at all 
contacts. The method to carry out this calculation is of- 
ten called the solver in this context, which is commonly 
an iterative scheme, as the one we describe in the follow- 
ing section. In order to make our points very clear, we 
perform an analytical investigation for a one dimensional 
example but can prove the existence of the discovered ef- 
fect in two dimensions as well. 
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A. The dynamical equations 

As CD is designed to obey excluded volume constraints 
exactly, particle collisions lead to discontinuous velocity 
changes ("shocks"), i.e. to nonsmooth mechanics [ p^ . 
Therefore, higher order terms than employed in the Euler 
integration scheme are of no use. For the ith particle's 
positions Xi, this reads 




ext 



FIG. 1: A pair of particles 



Xi{t + At) = Xi{t) -^v^{t + At) At , 
and correspondingly for its velocity Vi we have 



v,{t + At) = v,{t) + 



F^it + At) 



At 



(1) 



(2) 



where Fi is the total force acting on the particle, m its 
mass and At is the time step. 

A remark about the Euler scheme (|^) being implicit is 
in order: Whereas in conventional MD the choice of eval- 
uating the force for the previous or for the new configu- 
ration (i.e. at t or t + At respectively) is merely a matter 
of stability of the integration scheme (cf. e.g. [|9|), the 
constraints in CD can only be imposed on the new con- 
figuration. In this sense, the integration scheme of CD is 
inevitably of implicit type. (Taking into account the con- 
figurational change during a time step consistently {fully 
implicit integration |^) leads to difhculties which are 
analogous to implicite schemes in MD, when the forces 
have to be evaluated for the yet unknown new configura- 
tion. In one dimension, though, these difficulties do not 
arise.) 



1. One contact 

We now turn to the force on the I'th particle, Fi, occur- 
ring in Eq. (|^) : In order to determine it one has to know 
the contact forces between the grains. In CD they are 
calculated from the condition that the constraints must 
not be violated. In one dimension and if one disregards 
rotations, this is simply the excluded volume contraint. 
To give a specific example, let us consider two parti- 
cles with equal masses m subjected to constant external 
forces (Fig. |^). A contact force R (which is the reaction 
force due to the constraint) is active only, if interpenetra- 
tion needs to be prevented. Otherwise, i.e. if the gap g 
would remain non-negative (no overlap) anyway, it takes 
on its minimal value, R = (this is expressed in the CD 
literature as Signorini's condition). Without the repul- 
sion R between the particles the gap at the end of a time 
step would be given by 



g' = 9t + {v2,t - vi^)At 



^cxt 



At^ 



(3) 



according to (|I|) and (g). However, the excluded vol- 
ume constraint requires that gt+At — ma-xj^', 0}. In or- 
der that this results from (|^) and (^, the contact force 



Rt+At = max{i?', 0} must be taken into account, where 



R' 



mg 
2At^ 
m f-gt 



2At \ At 



+ Vit - V2t + 



(4) 
(5) 



Note that this scheme corresponds to a completely inelas- 
tic collision (i.e. with the so-called restitution coefficient 
being zero), which is accomplished in two time steps: At 
first the gap closes, in the next step also the relative ve- 
locity vanishes. (Finite restitution coefficients can also 
be incorporated into this algorithm po| .) 

The above determination of the contact force has a 
drawback, though: If gt < occurs due to a previous 
inaccuracy, then the elimination of this overlap is accom- 
panied by a surplus of kinetic energy. Therefore, mostly 
the quasi-inelastic shock formula p] 



2At 



pos 

-9t 
At 



+ Vif - V2,t + 



(6) 



is used instead of (||), where gf°^ = max{0,(?t}. That 
means, negative gaps are treated differently from Eq. (|^) 
in such a way that an already existing overlap is not elim- 
inated but only its further growth is inhibited. Hence the 
inelastic shock law (||) is in a way even "more inelastic" 
than the original law (^, because it avoids overlap cor- 
recting impulses which could destroy stable equilibrium 
states. 



2. Many contacts 

We now address the question how to solve the problem 
of the constraint forces if we consider not only one, but 
many contacts at the same time. Fig. |^ shows such a 
system, where Ri denotes the contact force between the 




FIG. 2: A multi-contact situation in an ID array. External 
forces act only on particles far away from those shown. Each 
particle is subjected only to the contact forces of the adjacent 
ones. 
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particles i and i + 1, and for the sake of simplicity we 
have no external force acting on them. Furthermore, we 
will concentrate on the situation where the neighboring 
pairs are permanently in contact (i.e. gf°^ = 0) which 
corresponds to compressed dense packings. 

For the ith contact in this setup, Ri-i and Ri+i play 
the role of the external forces in Eq. (||), though they 
are not constant, their values are not known for the next 
time step: 



Ri,t+At — 2^ {''^ht 



Ri-ld+At + Ri+l.t+At 



(7) 



Obviously the contact force is coupled to the neighboring 
contacts and through those to further contacts. In a sim- 
ilar way in higher dimensions, large numbers of contacts 
are coupled within clusters defined by the contact net- 
work. Hence, the determination of the proper reaction 
forces becomes a global problem. 

A standard way in CD is to solve this problem itera- 
tively: In one iteration step we calculate the forces ac- 
cording to the constraint conditions pretending that the 
corresponding neighboring contacts already exhibit the 
right forces. In that way the process traverses the list 
of contacts many times until satisfactory convergence is 
reached (for the question of the convergence cf. the works 

For our one-dimensional chain of particles, this means 
that Eq. (0) simply gets the meaning of an assignment 
of the right hand side to the left hand side, and one it- 
eration step consists of applying this assignment sequen- 
tially once to each contact. The order of this sequence is 
preferably random (with the pattern changing for every 
sweep) , in order not to create any bias in the information 
spreading With each sweep a globally consistent so- 
lution is approached, until finally a chosen convergence 
criterion is fulfilled. 

In the next section we shall show, what kind of conse- 
quences the local update scheme presented here has on 
large time and length scales. 



III. THE LARGE SCALE DESCRIPTION 

In order to analyze the coarse grained behavior of the 
microscopic equations derived in the last section we can 
also regard them as the discretized form of a continuum 
description, making a treatment in terms of partial dif- 
ferential equations (PDF) possible. In order to obtain 
the corresponding PDEs, we consider the particle index 
i as space variable x and replace the differences of con- 
secutive quantities by derivatives, the error term for the 
first and second order derivatives being of first and sec- 
ond order, respectively. E.g. Vt+At — vt ^ Atdtv and 
Ri+i + Ri-i — 2i?i — > (Pd^R, where d is the particle 
diameter. 



A. The relaxation of the contact forces 

While the continuum versions of the updates (|^) and 
(2) can be obtained straight forwardly, the force change 
(7) lacks a time variable, for during the force- iteration, 
being just a calculation, no physical time passes. Hence, 
to be able to describe this force-development as well, let 
us introduce a fictitious time t* with time interval At* 
for one iteration-step. With this, the continuum version 
of Eq. ^) reads: 

dt'R^ Dd^R- Pd^v (8) 
with D = (9) 



(3 



At* 
md 

^ AtAt* 



and 



1 



(10) 



(11) 



This analytic form clearly reveals the nature of the itera- 
tion loop: The reaction forces relax towards the solution 
in a diffusive way. (Note that the dxV term is constant 
in i*, it only depends on x.) 

The introduction of the constant q reflects a subtlety 
regarding the sequential character of the update dis- 
cussed in the appendix. In fact, since the PDF (^ de- 
scribes the change of the whole field R{x) at once, given 
its actual value at time t* , it corresponds to a parallel 
update (in the sense that the right hand side of Fq. (|^) 
always employs the values Ri from the beginning of the 
iteration sweep, not the freshly updated ones). In ap- 
pendix ^ we shall show, though, that a random sweep 
update instead of a parallel one only renormalizes the 
value of q to about 0.8 while leaving the form of the 
PDF untouched. 



B. Sound waves 

To connect the velocity update, whose continuum ver- 
sion follows immediately from Fq. (^ as 



dtV = dxR , 

TO 



(12) 



to the force update, we must relate the "iteration time" 
t* to the physical time t. Although, depending on the 
convergence criterion, there can be in principle a varying 
number of iterations during one physical time step At, 
we assume for simplicity this number Ni being fixed. 
(Actually, in practice this crude "criterion" is sometimes 
applied.) 

Hence, with At = NjAt*, we can express everything 
in terms of the physical time: 

dtR = Dd^R - PO^V , 



and D = qNi — 
At 

^ = ''''A^ 



(13) 
(14) 

(15) 
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With the equations (y_2|) and (|13|) we obtained two cou- 
pled PDEs. We can combine them to arrive at a wave 
equation with an additional damping term: 

djR = c^dlR + dt [DdlR] (16) 

The sound velocity appearing is of finite value: 

d 



VqNi 



At 



(17) 



This equation indicates that the CD simulation of 
the particle chain, as presented in the previous sec- 
tion, can lead to sound propagation like in an elastic 
medium, which however contradicts the conception of 
perfect rigidity. The constraint conditions applied at the 
contacts should in principle prohibit overlaps, i.e. pro- 
hibit elastic deformation of the grains. It can be seen 
that this deviation from the perfect rigidity enters at the 
force relaxation: A finite number of iterations means a 
finite range for the information spreading and thus yields 
systematic errors in the calculated reaction forces. As a 
consequence, the finite Nj involves soft particles and a 
finite sound velocity c ~ y^Nj. Note that in the limit of 
an infinite Nj, the exact value of the forces is reached, 
which corresponds to the case c — *■ oo, as it should be for 
rigid particles. 



1. Dispersion 

Performing a Fourier transformation on Eq. (p^, one 
obtains the properties of the different wave modes. The 
oscillation frequency uj of the wave number k is 



Lo (k) — k\ c 



D^k^ 



2. Numerical confirmation 

In order to confirm the results of this section, we per- 
formed the following numerical experiment: The start- 
ing configuration of the simulation consists of an array 
of 50 discs and an immobile wall, the geometry can be 
seen in Fig. ^. Initially the gap between the wall and 
the leftmost particle is one disc diameter (d), the gap 
between the particles is zero and the array has zero ve- 
locity. Starting from t = a constant external force 
[F°^^) is acting on the rightmost particle which acceler- 
ates the array towards the wall (only horizontal motion 
takes place). As simulation parameters we chose Nj — 40 
and F'''^' = 0.05 dmAt-^. 

The collision with the wall induces a relative motion of 
the grains and generates sound waves in the array. After 
a transient period the grains remain permanently in con- 
tact (the whole array is pressed against the wall by F°'^^). 
Since the different wave modes have different relaxation 
time, after a while only the largest wave length mode 
survives. This wave length is four times the system size 
because the wall represents a fixed boundary while the 
right side is free. Since the wave length is given, the oscil- 
lating frequency and the damping time can be calculated 
from Eq. ( p^ ) and Eq. (^^, respectively. For compari- 
son with the simulation we measured the motion of the 
rightmost particle. The expected motion is a damped 
oscillation 



x{t) = xo + ^exp {—t/r) sin [tot + i 



(21) 



(18) 



where the offset xq, the amplitude A and the phase shift 
(j) have to be fitted (in contrast to lu and r) for a compar- 
ison. In Fig. ^ the measured data (dots) and the fitted 
curve can be seen. It shows that the simulation is in good 
agreement with our continuum description. 



That means, lo (fc) becomes zero at a critical wave number 

1 



2c 

15 



(19) 



and waves with k larger than kc (short wave lengths) are 
over-damped. The damping time r(fc) for the oscillating 
modes is given by: 



(20) 



We derived the dispersion relation ( pq ) in the contin- 
uum limit which is a good approximation for small wave 
numbers, but not near to the border of the Brillouin zone 
[k-Qr — 27r/(i), where the effect of the spatial discreteness 
can be strong. However, increasing the number of the it- 
erations sufficiently, kc becomes small compared to k^^. 
Actually, for Ni > 10 the formula (|l^) works well not 
only for small wave numbers but for all the oscillating 
modes, as could be verified numerically. 



C. Global Elasticity 

It is instructive to compare our test-system to its sim- 
plest MD counterpart where the contact forces depend 
linearly on the local kinematic variables, i.e. the so called 
linear spring /dashpot model 

Ri = -K {xi+i - Xi- d) - J {vi+i - Vi) (22) 

with the spring stiffness k and the darnping coefficient 7. 
Employing again the updates (|^) and (||) for the positions 
and velocities, respectively, the continuum limit yields 



FIG. 3: The initial configuration of the numerical experiment 
for testing the properties of the sound waves. 
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FIG. 4: Damped oscillation in a Contact Dynamics simula- 
tion. The dots indicate the measured data: the position of the 
rightmost particle versus time (for details see the text). The 
line is an exponentially damped sinus function, where the fre- 
quency and the damping time is provided by our continuum 
model. 



the same type of PDE as Eq. ( |l^ ) with its coefficients 
being inherited from Eq. (22): 




dtR = —dlR - Kdd:,v 



(23) 



FIG. 5: The setup of a numerical experiment in two dimen- 
sions. A dense packing of 1000 discs is prepared in a container 
via compressing the system by means of the mobile upper 
wall. 



to time step), and therefore steps with a different "stiff- 
ness" are mixed during the integration of motion. Conse- 
quently, the behavior of the CD method is more complex 
in detail, but qualitatively the results for the constant Nj 
remain relevant also here. (E.g. the mechanism resulting 
in soft particles or the way how shock-waves can arise 
with finite velocity.) 



IV. 2D SIMULATION 



This allows us to relate the physical MD model parame- 
ters to the "technical" CD parameters: 



K = qm 



Ni 
At2 



and 



Ni 

7 = ?™^ 



(24) 



(25) 



This equivalence shows that on large scales the CD chain 
should behave identical to its MD counterpart, e.g. it 
will exhibit a global shrinkage proportional to an ex- 
ternal compressive load. Note that a real congruence 
can be expected only for the collective behavior but not 
on the level of the contacts. In the CD method, as ex- 
plained above, the contact forces are not related to the 
overlaps, which must merely be regarded as due to the 
incompleteness of the force-calculation and in fact are 
stochastic quantities because of our random update pro- 
cedure. Only on scales larger than the grain size, where 
the fluctuations of these local "deformations" are aver- 
aged out, the behavior can be smooth like in an elastic 
medium, as is shown in Fig. H. 



In sections [II B and [II C, our calculation was based 



on the assumption of a constant number of iterations 
for every time step, and due to this premise the analyti- 
cal treatment became simple and directly comparable to 
the corresponding simulation. We should keep in mind, 
though, that the application of a convergence criterion 
involves a fluctuating Nj (i.e. it can vary from time step 



After the analysis of the regular ID system, the im- 
portant question arises whether the behavior is similar 
in higher dimensions and less regular systems. Hence, 
we performed CD simulations with two-dimensional ran- 
dom packings of discs and observed the same "elastic" 
waves (even transversal modes were found) . 

The simulation presented here consists of 1000 discs 
with radii distributed uniformly between r^in and rmax — 
2?'minj the mass of each disc being proportional to its 
area. Fig. || shows the geometry: The base and the 
two side-walls are fixed while the upper piston is mo- 
bile. Starting from a loose state, we compressed the sys- 
tem and waited until the packing reached an equilibrium 
state (the compression force F applied on the piston was 
kept constant). The simulation was carried out without 
gravity and with a Coulomb friction coefficient of 0.05 
for all the disc-disc and disc- wall contacts (cf. [p3[). 

After the packing was relaxed completely, we generated 
sound waves by increasing the compression force abruptly 
to F + AF. After a transient period only one stand- 
ing wave mode survives (both the wavenumber vector 
and the collective motion are vertical), where the piston, 
representing a free boundary, oscillates with a relatively 
large amplitude. We measured the vertical position of 
the piston versus time and found that the data can be 
fitted by an exponentially damped sine function (Fig. |^). 
Here, in contrast to the ID case, also lu and r are fit 
parameters, since, due to the different geometry, the val- 
ues ( |l8| ) and ( |2C| ) cannot be adopted, but, because of the 
random structure of the system, a more complex treat- 
ment is required for a quantitative description. However, 
we checked the most important relation, namely that the 
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FIG. 6: Oscillations in a 2D simulation are similar to the ID 
case. Here the sound waves are generated in a random dense 
packing of discs. The dots are the measured position of the 
upper wall versus time (see Fig. ^ , while the curve is a fitted 
exponentially damped sinus function. 

scaling properties of lo and t remain valid also for the 2D 
random system, that is w ^ \/ Ni and r ~ NJ^ , which 
means that the artificial visco-elasticity of the particles 
depends on the number of the iterations in the same man- 
ner as we showed for the ID chain. 



number of particles n can be determined: One step of 
the iteration consists of as many force-updates as there 
are contacts, which is proportional to n. Therefore, the 
computational effort of one time step scales with the par- 
ticle number like nNj^ which is ^ "n? in 2D or ~ rv'/^ in 
3D. Therefore a large CD-simulation is computationally 
more costly than MD, where the computational effort 
scales like n. This is the price for simulating rigid parti- 
cles without getting elasticity artifacts, which cannot be 
done with MD. 

To avoid this super-linear scaling when dealing with 
large systems, we can also accept the finite stiffness by 
keeping Nj constant independently of n. Then, besides 
gaining a running time of order n, of course, elastic de- 
formations and sound waves can arise with an increasing 
number of the particles, and consequently they have to be 
monitored. We want to mention the idea, though, that 
in certain situations advantage can be taken of the arti- 
fact. E.g. when being applied deliberately, Coulombian 
friction can be combined with global elasticity easily; this 
way considerable computational time could be saved and 
even better performance than MD could be achieved. 
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V. DISCUSSION 

The artificial elasticity found in CD simulations was 
analyzed. We showed that the systematic errors of the 
force calculation can lead to a collective elastic behavior, 
even though single contacts are assumed perfectly rigid 
and perfectly inelastic. For the ID chain of particles, we 
could, starting from the "microscopic" laws, reproduce 
the numerical results analytically, including the depen- 
dence of the effective stiffness and viscous dissipation of 
the contacts on the computational parameters {Nj^At). 

Besides elucidating the origin of the elastic behaviour, 
the coarse grained description reveals important charac- 
teristics of the CD method, which were less obvious on 
the discrete level. We saw that using the iterative solver, 
the proper contact forces are approached in a diffusion 
like manner which is a crucial information concerning the 
computational time. The conception of perfectly rigid 
particles requires that the calculated forces are consis- 
tent even for contacts f ar fro m each other. Therefore, 
the "diffusion length" V DAt must be larger than the 
linear system size i, which defines a lower boundary for 
the number of iterations of the order (L/c?)^. The same 
condition is obtained if we want to avoid other conse- 
quences of the effective softness. E.g. if we want the 
sound to travel a larger distance than L during one time 
step (i.e. c > L/At) or if we would hke all possible wave 
modes to be overdamped (i.e. with Ac > L no oscilla- 
tions). All these cases are equivalent, one is forced to 
apply a relatively large number of iterations: Nj ~ L^. 
Going further on this line, the scaling with respect to the 
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APPENDIX A: 



FOR THE RANDOM SWEEP 
UPDATE 



To find the proper value of the constant q appearing in 
two more things have to be taken into account: 



II] 



sec 

Firstly the sequential type of the force update, and sec- 
ondly that the order is random. The latter results in a 
stochastic force relaxation, i.e. the change of Ri in one it- 
eration sweep is a stochastic variable. In order to obtain 
a similar equation as (|^) we shall determine the average 
value (ARi), where the average for site i is meant to be 
taken over all possible update sequences. 

Before going any further, let us introduce a few nota- 
tions: 

• For N being the total number of contacts, the map- 
ping u : {1, . . . , N} — > {1, . . . , N} denotes the or- 
der of the update sequence; i.e. if the contact la- 
beled i is updated before j, then Ui < Uj. 

• Throughout this appendix, the notation Ri means 
the value from the beginning of the iteration sweep, 
the value at the end is Ri + ARi . 

• We define 5Ri as the change according to a parallel 
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update (cf. Eq. (0)), i.e. 



Ri-i + — 2Ri 



as opposed to the total change A_Ri 



(Al) 



Given that, it can easily be seen how Ai?^ depends on 
the update order. If e.g. the site i is updated earher than 
its neighbors (i.e. Ui-i > Ui < Wi+i), then ARi = SRi, 
but in the case Ui-i > m > Ui+i < we get 

Ai?i = 5Ri + SRi+i/2 (because for contact i, Eq. (Q) 
employs the already updated force at contact i + Sim- 
ilarly (but with less probability) even very far contacts 
can contribute to ARi, which can be summarized in the 
following way: 



{AR, 



5R^ 



{PrSRi-r +PrSRi+r) 



(A2) 



Here, is the probability that ARi contains informa- 
tion from the update of a contact at distance r, that is 
Pr — P{ui > Ui+i > . . . > Mi+r)- (This definition holds 
true for contributions from contacts with labels higher 
than i, but due to left-right symmetry the same value is 
inevitably obtained for the corresponding lower ones.) 

The value of pr can be obtained from the following 
combinatorial consideration: Given an index i and a dis- 
tance r, we can classify the set of all update orders into 



groups such that the sequences in one group differ only in 
the permutations of the elements Uj, i < j < i -\- r. Such 
a group contains {r + l)l sequences, but only one of them 
satisfies the condition Ui > u^+i > Mi+2 > . . . > Ui+r- 
Since all update sequences are equally probable, the value 
of Pr is equal to l/(r + 1)!. 

The factor {2^{r + 1)!)~^, relating the contacts i and 
i + r, decays faster than exponentially; already for r — 8, 
it drops below 10^®. Therefore, the sum in Eq. (A2) 
reaches only the immediate vicinity of contact i, such 
that, for our large wavelength considerations, the approx- 
imation SRi+r ~ SRi can be applied. This allows us to 
calculate the average change of the contact force: 

= (4V^-5) . (A3) 



Thus, it is shown that the random sweep results in a 
larger change of Ri than the parallel update. Eq. ( [A^ ) 
provides also the sought value of the parameter q as 



, = « 0.797 



(A4) 



which completes the continuum description given in 
sec. [IL 
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